##########################################################
# John Henderson
# Replication Data for: "Hookworm Eradication as a Natural 
#   Experiment for Schooling and Voting in the American 
#   South", Forthcoming in Political Behavior, May 20, 2017
# 
##########################################################
#
#  figure 4a.R
#  -- file produces exclusion sensitive plots for the RSC campaign
#
##########################################################

rm(list=ls())

library(Matching)
library(AER)                           
library(foreign)

setwd('~/Dropbox/hookworm/replication')
source('funs/lower.bound.R')
            
load('data/hookwormData.Rdata') 

##########################################################  
# data required 
#  - mMat :: matching matrix used in matching and balancing
#  - oMat :: outcome matrix of county education and turnout rates
#  - zMat :: intervention matrix 
#  - pMat :: placebo matrix
# hookSouth :: combines all objects into one dataset, w/ additional covariates  
##########################################################

##########################################################
##########################################################
# BEFORE MATCHING
##########################################################
##########################################################

d1=oMat$school1920
d0=oMat$school1910    

y0=rowMeans(cbind(oMat$vote1916,oMat$vote1920))
y1=rowMeans(cbind(oMat$vote1928,oMat$vote1932))

tr=zMat$z_campaign

d=d1-d0
y=y1-y0
d_dd=c(d0,d1)
y_dd=c(y0,y1)
tr_dd=c(tr,tr) 

xmat=mMat[,-c(1:3)][,1:19]
xmat=as.matrix(rbind(xmat,xmat))

year=c(rep(0,length(tr)),rep(1,length(tr)))
temps=as.data.frame(cbind(d_dd,y_dd,tr_dd,year))   

##########################################################
##########################################################
# AFTER MATCHING
##########################################################
##########################################################

mMat=bMat=mMat[,-c(1:3)][,1:19]

# pre-sets for 'lower.bound' fitfunc :: lower.bound only improves if total balance improves
floors=NULL					# :: allows to preset values at the floor overall
floor.values=NULL           # :: allows to preset values at the floor for each covariate
reweighted = -1 			# :: fits over a single summary measure of balance
exacts=array(0,ncol(mMat))  # :: exact matching for covariates 
exacts[1]=0    

starts=c(0.8639389,12.5338719,14.1059201,232.6227015,13.3530437,0.8680492,809.5765534,2.3393068,1.8030039,0.9325587,44.5785061,1.5006086,65.5576631,2.8099156,9.4050025,43.4617957,453.9672492,357.1535636,388.9175985)     

set.seed(1005)
genout=GenMatch(starting.values=starts,exact=exacts,MemoryMatrix=T,fit.func=lower.bound,Tr=tr,X=mMat,BalanceMatrix=bMat,M=1,pop.size=2,max.generations=5,wait.generations=3,hard.generation.limit=T,estimand="ATT",ties=T)

mout=(Match(Tr=tr,X=mMat,estimand="ATT",M=1,Weight.matrix=genout,exact=exacts,ties=T))

trs=c(rep(1,length(mout$index.treated)),rep(0,length(mout$index.control)))
xcovar=c(mout$index.treated,mout$index.treated)

dtr=d[mout$index.treated]
dct=d[mout$index.control]
ytr=y[mout$index.treated]
yct=y[mout$index.control]
w=mout$weights
trt=tr[mout$index.treated]
tct=tr[mout$index.control]

ymat=c(ytr,yct)
dmat=c(dtr,dct)
zmat=c(trt,tct)
wmat=c(w,w)
      
xb=ivreg(dmat~zmat,weights=wmat)$coef[2] 
ivb=(ivreg(ymat~dmat | zmat,weights=wmat))$coef[2]

gammas=seq(from=0,to=ivb*xb,by=.00025)
cfs=ci1=ci2=array(0,length(gammas))
for(i in 1:length(gammas)){
	yadj1 = ymat - gammas[i]*zmat
	yadj2 = ymat + gammas[i]*zmat	
	sum1=summary(ivreg(yadj1~dmat | zmat,weights=wmat)) 	
	sum2=summary(ivreg(yadj2~dmat | zmat,weights=wmat))

   ci1[i]=sum1$coef[2,1]-1.96*sum1$coef[2,2]   
   ci2[i]=sum2$coef[2,1]+1.96*sum2$coef[2,2]	
}

                       
gtp=((gammas/xb)/ivb)[min(which(ci1<0))]   
                

pdf('~/Dropbox/hookworm/replication/figures/figure_4a_campaignExcluded.pdf')  
par(cex=1.2)   
plot(x=-100,y=-100,xlim=c(0,1),ylim=c(-.75,3.5),ylab='',xlab='')
labNames=c(' IV Coefficient',' Ratio of Coefficients')
mtext(2,line=2.75,text=substitute(paste(beta,'(',gamma,')',nn), list(nn=labNames[1])) ,cex=1.3)   
mtext(1,line=2.75,text=substitute(paste(gamma/hat(beta),nn), list(nn=labNames[2])) ,cex=1.3)    
lines((gammas/xb)/ivb,ci1,col='red',lty=2,lwd=3)
lines((gammas/xb)/ivb,ci2,col='red',lty=2,lwd=3)  
abline(h=0,lty=2,col='grey60')
lines(x=c(gtp+.02,gtp+.135),y=c(0.045,.75))
text(x=gtp+.145,y=.9,label=round(gtp,digits=2),cex=1.2) 
text(x=.135,y=3.2,label=(expression(alpha[CI]==0.95)),cex=1.2)      
dev.off()
          

pdf('~/Dropbox/hookworm/replication/figures/figure_4a_campaignExcluded_bw.pdf')  
par(cex=1.2)   
plot(x=-100,y=-100,xlim=c(0,1),ylim=c(-.75,3.5),ylab='',xlab='')
labNames=c(' IV Coefficient',' Ratio of Coefficients')
mtext(2,line=2.75,text=substitute(paste(beta,'(',gamma,')',nn), list(nn=labNames[1])) ,cex=1.3)   
mtext(1,line=2.75,text=substitute(paste(gamma/hat(beta),nn), list(nn=labNames[2])) ,cex=1.3)    
lines((gammas/xb)/ivb,ci1,col='black',lty=2,lwd=3)
lines((gammas/xb)/ivb,ci2,col='black',lty=2,lwd=3)  
abline(h=0,lty=2,col='grey60')
lines(x=c(gtp+.02,gtp+.135),y=c(0.045,.75))
text(x=gtp+.145,y=.9,label=round(gtp,digits=2),cex=1.2) 
text(x=.135,y=3.2,label=(expression(alpha[CI]==0.95)),cex=1.2)      
dev.off()

#END